% Script for the structural analysis of BFW data in 
% Marra, Marleen (2024), "Estimating an Auction Platform Game with Two-sided Entry",
% Journal of Political Economy, (in print July 2024)
% 
clc
clear
rng('default');
matlabversion = version
date = date


% ==================== Main specification in paper ============================= %
% While some "options" below can be changed, the code is only guaranteed for the selected values. 

% FIRST STAGE ----------------------------------------- %
% Check "firststage_modelselection.m" for why these variables are selected
PAR.Firststage_model = ['log_highestbid_bot~specialformat+onebottle+numberbottles+numberbottles2+caseof6+caseof12+'...
    'inwarehouse+'...
    'onlycollect+returnsaccepted+insuranceincl+'...
    'filllevel_ordered+winetype+'...
    'popularregions+'...
    'deliverstouk+canbecollected+'...
    'paymentbank+paymentpaypal+paymentcheque+paymentcash+'...
    'shiproyalmail+shipparcelforce+shipfast+monthended_group+'...
    'dutyestimate+vatestimate+ship+hasratings+sellerratings+sellerratings2+'...
    'desc_primeur+desc_delivery+desc_parker+desc_numberwords'];

% PARAMETERS for estimation --------------------------- %
% Some empirical choices
PAR.trimQhigh   = 0.02; % How much of estimated first-stage quality trim off top and bottom (drop observations)
PAR.truncu0high = 0.02; % How much of implied seller value to trim off bottom and top for estimation (top gets replaced with u0star in iterations >1)
PAR.truncu0main = 0.01; % How much of implied seller value to trim off bottom and top for estimation (top gets replaced with u0star in iterations >1) 
PAR.trimQmain   = 0.01; % How much of estimated first-stage quality trim off top and bottom (drop observations)

% Options for simulation of entry equilibrium
PAR.maxn              = 15;             % Reasonable upper bound on number of bidders per listing
PAR.sims              = 250;            % To calculate HV, SHV, then interpolate on fine grid 1000*sims, then 100*sims MC draws from values exceeding rhat
PAR.v0probsfine       = unique(sort(random('unif',0.0001,0.9999,[500,1])));  % below Determines grid of v0star - used to calculate bidder surplus, lambdastar, v0star
PAR.v0probs           = unique([0.0001;sort(randsample(PAR.v0probsfine,48));0.9999]); % then interpolate on finer grid
PAR.opts              = optimset('MaxFunEvals',200,'MaxIter',100); %fmincon in Estimation_nestedfp.m
PAR.opts.Display      ='off';
PAR.optsv0star        = optimset('Display','off','TolX',1e-6,'TolFun',1e-6,'UseParallel',false); 
PAR.optslambdastar    = optimset('TolX',1e-6,'TolFun',1e-6,'UseParallel',false);
PAR.probsMC           = unique(sort([0.01;random('unif',0,1,[PAR.sims,1]);0.99]));
PAR.MC                = sort(random('unif',min(PAR.probsMC),max(PAR.probsMC),100*PAR.sims,1));

% NPL convergence
PAR.NPLtol            = 0.01/4; %convergence of max(abs([ F_V0(v0*) , eSo, thetas ]))
PAR.NPLmaxit          = 100;

% Parameterization DISTRIBUTION private values bidders & sellers --- %
PAR.dist              = 'normal';       % Parameterization idiosyncratic value distributions
PAR.generalized       = 1;              % Use generalized Gaussian
PAR.th0_main          = [4,0.5,-0.1];   % Starting parameters for main sample
PAR.th0_high          = [5,0.5,-0.1];   % Starting parameters for high-value sample
PAR.lb                = [-Inf,0,-2.5];  % The last value for shape parameter if generalized=1
PAR.ub                = [Inf,Inf,2.5];  % The last value for shape parameter ifgeneralized=1
PAR.add0topois        = 1;              % USe generalized Poisson distribution

% Bootstrapping
PAR.nrbootstraps      = 250;            % Number of bootstrap iterations
PAR.thres_mus         = 50;             % Threshold of seller mu_s above which estimates are considered "blown up" (and excluded) based on unfortunate sample draw
PAR.values_gridpoints = 100;            % Number of gridpoints for CDF plots estimated value distributions
PAR.bootstrap         = 0;              % Automatically changes to 1 when doing bootstraps, just so that model fit is based on main estimates
dobootstrap           = 1;              % Set to 1 if want to compute bootstrap CI as well

% Model fit
PAR.v0probs_fit       = random('unif',0.0001,0.9999,[5000,1]);        % large nr draws - also take expectations over draws of Quality     
PAR.vprobs_fit        = random('unif',0.0001,0.9999,[1000,PAR.maxn]); % large nr draws - also take expectations over draws of Quality  


% ======================== Main specification in paper  ================================== %

% Load data
data                  = readtable('BFWauctiondata.csv','TreatAsEmpty','NA');

%% Estimation 
% main sample (0-200]
PAR.highvalue = 0; 
[Results,ESTIntermediates,ESTPars,ESTValues,ESTNetworkEffects,ESTNPLiterations] = estimation_func(data,PAR); % point estimates
[BSResults,BSidx,noblow,BS,BSResultsorig] = bootstrap_func(data,PAR,Results,'../Output/'); % bootstrapped standard errors
Modelfit = modelfit_func(PAR,Results,'../Output/'); % Model fit checks 

%  high-value sample (200-800)
PAR.highvalue = 1;
[Results_high,ESTIntermediates_high,ESTPars_high,ESTValues_high,ESTNetworkEffects_high,ESTNPLiterations_high]    = estimation_func(data,PAR); % point estimates
[BSResults_high,BSidx_high,noblow_high,BS_high,BSResultsorig_high] = bootstrap_func(data,PAR,Results_high,'../Output/'); % bootstrapped standard errors
Modelfit_high = modelfit_func(PAR,Results_high,'../Output/');% Model fit checks 

% Produce Table 5 - Estimation results (labels of columns in 3/4th lines, and brackets to be added manually)
RES=round([Results.mub,Results.kb,Results.shapeb,Results.mus,Results.ks,Results.shapes;...
                   BS.BS_SE(9),BS.BS_SE(10),BS.BS_SE(13),BS.BS_SE(11),BS.BS_SE(12),BS.BS_SE(14);...
                   Results.eBo_posr,Results.eBo_nor,Results.eSo,Results.sharenobidders,Results.screeningvalue,NaN; ...
                   BS.BS_SE(6),BS.BS_SE(7),BS.BS_SE(8),BS.BS_SE(21),BS.BS_SE(15),NaN],3);...
        matrix2latex(RES, '../Output/Table5_PanelI.tex',...
            'ColumnLabels',[{'mub', 'varb','shapeb','mus','vars','shapes'} ],...
            'RowLabels',[{'Values est'},{'Values se'},{'Entry est'},{'Entry se'}])
 
RESH=round([Results_high.mub,Results_high.kb,Results_high.shapeb,Results_high.mus,Results_high.ks,Results_high.shapes;...
                   BS_high.BS_SE(9),BS_high.BS_SE(10),BS_high.BS_SE(13),BS_high.BS_SE(11),BS_high.BS_SE(12),BS_high.BS_SE(14);...
                   Results_high.eBo_posr,Results_high.eBo_nor,Results_high.eSo,Results_high.sharenobidders,Results_high.screeningvalue,NaN; ...
                   BS_high.BS_SE(6),BS_high.BS_SE(7),BS_high.BS_SE(8),BS_high.BS_SE(21),BS_high.BS_SE(15),NaN],3);...
        matrix2latex(RESH, '../Output/Table5_PanelII.tex',...
            'ColumnLabels',[{'mub', 'varb','shapeb','mus','vars','shapes'} ],...
            'RowLabels',[{'Values est'},{'Values se'},{'Entry est'},{'Entry se'}])
                       

%% Policy simulations
% Baseline simulations of platform model at estimated parameters
shut_bidderentry = 0;
shut_sellerentry = 0;
PAR.highvalue    = 0; 
CFBaseline       = Simulate_CF(Results,PAR,Results,shut_bidderentry,shut_sellerentry);
PAR.highvalue    = 1; 
CFBaseline_high  = Simulate_CF(Results_high,PAR,Results_high,shut_bidderentry,shut_sellerentry);

% Counterfactuals 
[CF1,CF4,CF5,CF6,CF7,CF10] = counterfactual_func(Results,Results_high,CFBaseline,CFBaseline_high,PAR,'../Output/',1);%last entry = 0 to only produce figures based on already executed counterfactual simulations


save('../Output/All_results')








 

